//
// Created by Ivor_ on 2022/2/25.
//

#include "monteCarlo.h"

double sigmoid(double xx) {
    return 1./(1. + exp(-xx));
}

StdSqrNet createNetworks_stdSqr(int scales) {
    srand(time(NULL));
    StdSqrNet ssn;
    ssn.scales = scales;
    (ssn.network) = (Node*)calloc(scales*scales, sizeof(Node));
    for (int i = 0; i < scales; ++i) {
        for (int j = 0; j < scales; ++j) {
            if (rand()%2) {
                (ssn.network)[i*scales + j].attribute = 'c';
            } else  {
                (ssn.network)[i*scales + j].attribute = 'd';
            }
            (ssn.network)[i*scales + j].neighbors = 4;
            (ssn.network)[i*scales + j].neighbor = (Node**)calloc(4, sizeof(Node*));
            (ssn.network)[i*scales + j].neighbor[0] = &((ssn.network)[((i+scales-1)%scales)*scales + j]);
            (ssn.network)[i*scales + j].neighbor[1] = &((ssn.network)[((i+1)%scales)*scales + j]);
            (ssn.network)[i*scales + j].neighbor[2] = &((ssn.network)[i*scales + (j+scales-1)%scales]);
            (ssn.network)[i*scales + j].neighbor[3] = &((ssn.network)[i*scales + (j+1)%scales]);
        }
    }
    return ssn;
}

StdTreeNet createNetworks_stdTree(int layers) {
    srand(time(NULL));
    StdTreeNet stn;
    stn.layers = layers;
    (stn.network) = (Node*)calloc((size_t)(pow(3., layers+1.) - 1.), sizeof(Node));
    for (int i = 0; i <= layers; ++i) {
        for (int j = 0; j < (int)(2.*pow(3., i)); ++j) {
            if (rand()%2) {
                (stn.network)[(int)pow(3., i) - 1 + j].attribute = 'c';
            } else  {
                (stn.network)[(int)pow(3., i) - 1 + j].attribute = 'd';
            }
            if (i < layers) {
                (stn.network)[(int)pow(3., i) - 1 + j].neighbors = 4;
                (stn.network)[(int)pow(3., i) - 1 + j].neighbor = (Node**)calloc(4, sizeof(Node*));
                if (i == 0) {
                    (stn.network)[(int)pow(3., i) - 1 + j].neighbor[0] = &((stn.network)[1-j]);
                } else {
                    (stn.network)[(int)pow(3., i) - 1 + j].neighbor[0] = &((stn.network)[(int)pow(3., i-1.) - 1 + j/3]);
                }
                (stn.network)[(int)pow(3., i) - 1 + j].neighbor[1] = &((stn.network)[(int)pow(3., i+1.) - 1 + 3*j + 0]);
                (stn.network)[(int)pow(3., i) - 1 + j].neighbor[2] = &((stn.network)[(int)pow(3., i+1.) - 1 + 3*j + 1]);
                (stn.network)[(int)pow(3., i) - 1 + j].neighbor[3] = &((stn.network)[(int)pow(3., i+1.) - 1 + 3*j + 2]);
            } else {
                (stn.network)[(int)pow(3., i) - 1 + j].neighbors = 1;
                (stn.network)[(int)pow(3., i) - 1 + j].neighbor = (Node**)calloc(1, sizeof(Node*));
                (stn.network)[(int)pow(3., i) - 1 + j].neighbor[0] = &((stn.network)[(int)pow(3., i-1.) - 1 + j/3]);
            }
        }
    }
    return stn;
}

ArbGraph createNetworks_fourDegreeRing(int nodes, int gap) {
    srand(time(NULL));
    ArbGraph fdr;
    fdr.nodes = nodes;
    (fdr.network) = (Node*)calloc(nodes, sizeof(Node));
    for (int i = 0; i < nodes; ++i) {
        if (rand()%2) {
            fdr.network[i].attribute = 'c';
        } else {
            fdr.network[i].attribute = 'd';
        }
        fdr.network[i].neighbors = 4;
        fdr.network[i].neighbor = (Node**)calloc(4, sizeof(Node*));
        fdr.network[i].neighbor[0] = &(fdr.network[(i+1+nodes)%nodes]);
        fdr.network[i].neighbor[1] = &(fdr.network[(i-1+nodes)%nodes]);
        fdr.network[i].neighbor[2] = &(fdr.network[(i+gap+nodes)%nodes]);
        fdr.network[i].neighbor[3] = &(fdr.network[(i-gap+nodes)%nodes]);
    }
    return fdr;
}

ArbGraph ssnToAg(StdSqrNet ssn) {
    ArbGraph ag;
    ag.nodes = ssn.scales*ssn.scales;
    ag.network = ssn.network;
    return ag;
}

ArbGraph stnToAg(StdTreeNet stn) {
    ArbGraph ag;
    ag.nodes = (int)(pow(3., stn.layers) - 1.);
    ag.network = stn.network;
    return ag;
}

void randRelinkGraph(ArbGraph* pag, double possibility) {
    srand(time(NULL));
    for (int i = 0; i < (*pag).nodes; ++i) {
        for (int j = 0; j < (*pag).network[i].neighbors; ++j) {
            if (1.*rand()/RAND_MAX < possibility) {
                int objLoc;
                do {
                    objLoc = (abs(rand()) + RAND_MAX*abs(rand()))%(*pag).nodes;
                } while (objLoc != i);
                Node** tempNodeArray = (Node**)calloc(((*pag).network[i].neighbor[j])->neighbors - 1, sizeof(Node*));
                for (int k = 0, l = 0; k < ((*pag).network[i].neighbor[j])->neighbors; ++k) {
                    if (((*pag).network[i].neighbor[j])->neighbor[k] != &((*pag).network[i])) {
                        tempNodeArray[l] = ((*pag).network[i].neighbor[j])->neighbor[k];
                        l++;
                    }
                }
                free(((*pag).network[i].neighbor[j])->neighbor);
                ((*pag).network[i].neighbor[j])->neighbors -= 1;
                ((*pag).network[i].neighbor[j])->neighbor = (Node**)calloc(((*pag).network[i].neighbor[j])->neighbors, sizeof(Node*));
                for (int k = 0; k < ((*pag).network[i].neighbor[j])->neighbors; ++k) {
                    ((*pag).network[i].neighbor[j])->neighbor[k] = tempNodeArray[k];
                }
                free(tempNodeArray);
                (*pag).network[i].neighbor[j] = &((*pag).network[objLoc]);
                tempNodeArray = (Node**)calloc((*pag).network[objLoc].neighbors, sizeof(Node*));
                for (int k = 0; k < (*pag).network[objLoc].neighbors; ++k) {
                    tempNodeArray[k] = (*pag).network[objLoc].neighbor[k];
                }
                free((*pag).network[objLoc].neighbor);
                (*pag).network[objLoc].neighbors += 1;
                (*pag).network[objLoc].neighbor = (Node**)calloc((*pag).network[objLoc].neighbors, sizeof(Node*));
                for (int k = 0; k < (*pag).network[objLoc].neighbors - 1; ++k) {
                    (*pag).network[objLoc].neighbor[k] = tempNodeArray[k];
                }
                free(tempNodeArray);
                (*pag).network[objLoc].neighbor[(*pag).network[objLoc].neighbors - 1] = &((*pag).network[i]);
            }
        }
    }
}

void freeSSN(StdSqrNet ssn) {
    for (int i = 0; i < ssn.scales*ssn.scales; ++i) {
        free(ssn.network[i].neighbor);
    }
    free(ssn.network);
}

void freeSTN(StdTreeNet stn) {
    for (int i = 0; i < (int)(pow(3., stn.layers) - 1.); ++i) {
        free(stn.network[i].neighbor);
    }
    free(stn.network);
}

void freeAG(ArbGraph ag) {
    for (int i = 0; i < ag.nodes; ++i) {
        free(ag.network[i].neighbor);
    }
    free(ag.network);
}

double ssnMcApproximate(StdSqrNet* pStdSqrNet, double cost) {
    srand(time(NULL));
    double receive = cost + 1.;
    double density[2] = {0., 0.};
    for (int i = 0; i < (*pStdSqrNet).scales*(*pStdSqrNet).scales; ++i) {
        density[1] += 1.*(((*pStdSqrNet).network)[i].attribute == 'c');
    }
    density[1] /= pow(1.*(*pStdSqrNet).scales, 2.);
    int iterCount = 0;
    do {
        for (int i = 0; i < (*pStdSqrNet).scales*(*pStdSqrNet).scales; ++i) {
            int randObj[2] = {abs((rand() + RAND_MAX*rand())%((*pStdSqrNet).scales*(*pStdSqrNet).scales)), 0};
            randObj[1] = abs(rand()%(*pStdSqrNet).network[randObj[0]].neighbors);
            if ((*pStdSqrNet).network[randObj[0]].attribute != (*pStdSqrNet).network[randObj[0]].neighbor[randObj[1]]->attribute) {
                double go, gn;
                go = - cost * ((*pStdSqrNet).network[randObj[0]].attribute == 'c') * (*pStdSqrNet).network[randObj[0]].neighbors;
                for (int j = 0; j < (*pStdSqrNet).network[randObj[0]].neighbors; ++j) {
                    go += receive*((*pStdSqrNet).network[randObj[0]].neighbor[j]->attribute == 'c');
                }
                gn = - cost * ((*pStdSqrNet).network[randObj[0]].neighbor[randObj[1]]->attribute == 'c') * (*pStdSqrNet).network[randObj[0]].neighbor[randObj[1]]->neighbors;
                for (int j = 0; j < (*pStdSqrNet).network[randObj[0]].neighbor[randObj[1]]->neighbors; ++j) {
                    gn += receive*((*pStdSqrNet).network[randObj[0]].neighbor[randObj[1]]->neighbor[j]->attribute == 'c');
                }
                double possibility = sigmoid((gn - go)/.1);
                if (1.*abs(rand())/RAND_MAX < possibility) {
                    (*pStdSqrNet).network[randObj[0]].attribute = (*pStdSqrNet).network[randObj[0]].neighbor[randObj[1]]->attribute;
                }
            }
        }
        density[0] = density[1];
        for (int i = 0; i < (*pStdSqrNet).scales*(*pStdSqrNet).scales; ++i) {
            density[1] += 1.*(((*pStdSqrNet).network)[i].attribute == 'c');
        }
        density[1] /= pow(1.*(*pStdSqrNet).scales, 2.);
        iterCount++;
    } while (fabs(density[0] - density[1]) > 1.e-6 || iterCount < (*pStdSqrNet).scales);
    return (density[1]);
}

double stnMcApproximate(StdTreeNet* pStdTreeNet, double cost) {
    srand(time(NULL));
    double receive = cost + 1.;
    double density[2] = {0., 0.};
    for (int i = 0; i < (int)(pow(3., (*pStdTreeNet).layers+1.) - 1.); ++i) {
        density[1] += 1.*(((*pStdTreeNet).network)[i].attribute == 'c');
    }
    density[1] /= (pow(3., (*pStdTreeNet).layers+1.) - 1.);
    int iterCount = 0;
    do {
        for (int i = 0; i < (int)(pow(3., (*pStdTreeNet).layers+1.) - 1.); ++i) {
            int randObj[2] = {abs((rand() + RAND_MAX*rand())%(int)(pow(3., (*pStdTreeNet).layers+1.) - 1.)), 0};
            randObj[1] = abs(rand()%((*pStdTreeNet).network[randObj[0]].neighbors));
            if ((*pStdTreeNet).network[randObj[0]].attribute != (*pStdTreeNet).network[randObj[0]].neighbor[randObj[1]]->attribute) {
                double go, gn;
                go = - cost * ((*pStdTreeNet).network[randObj[0]].attribute == 'c') * (*pStdTreeNet).network[randObj[0]].neighbors;
                for (int j = 0; j < (*pStdTreeNet).network[randObj[0]].neighbors; ++j) {
                    go += receive*((*pStdTreeNet).network[randObj[0]].neighbor[j]->attribute == 'c');
                }
                gn = - cost * ((*pStdTreeNet).network[randObj[0]].neighbor[randObj[1]]->attribute == 'c') * (*pStdTreeNet).network[randObj[0]].neighbor[randObj[1]]->neighbors;
                for (int j = 0; j < (*pStdTreeNet).network[randObj[0]].neighbor[randObj[1]]->neighbors; ++j) {
                    gn += receive*((*pStdTreeNet).network[randObj[0]].neighbor[randObj[1]]->neighbor[j]->attribute == 'c');
                }
                double possibility = sigmoid((gn - go)/.1);
                if (1.*abs(rand())/RAND_MAX < possibility) {
                    (*pStdTreeNet).network[randObj[0]].attribute = (*pStdTreeNet).network[randObj[0]].neighbor[randObj[1]]->attribute;
                }
            }
        }
        density[0] = density[1];
        for (int i = 0; i < (int)(pow(3., (*pStdTreeNet).layers+1.) - 1.); ++i) {
            density[1] += 1.*(((*pStdTreeNet).network)[i].attribute == 'c');
        }
        density[1] /= (pow(3., (*pStdTreeNet).layers+1.) - 1.);
        iterCount++;
    } while (fabs(density[0] - density[1]) > 1.e-6 || iterCount < (int)(pow(3., (*pStdTreeNet).layers+1.) - 1.));
    return (density[1]);
}

double agMcApproximate(ArbGraph* pArbGraph, double cost) {
    srand(time(NULL));
    double receive = cost + 1.;
    double density[2] = {0., 0.};
    for (int i = 0; i < (*pArbGraph).nodes; ++i) {
        density[1] += 1.*(((*pArbGraph).network)[i].attribute == 'c');
    }
    density[1] /= 1.*(*pArbGraph).nodes;
    int iterCount = 0;
    do {
        for (int i = 0; i < (*pArbGraph).nodes; ++i) {
            int randObj[2] = {abs((rand() + RAND_MAX*rand())%(*pArbGraph).nodes), 0};
            randObj[1] = abs(rand()%(*pArbGraph).network[randObj[0]].neighbors);
            if ((*pArbGraph).network[randObj[0]].attribute != (*pArbGraph).network[randObj[0]].neighbor[randObj[1]]->attribute) {
                double go, gn;
                go = - cost * ((*pArbGraph).network[randObj[0]].attribute == 'c') * (*pArbGraph).network[randObj[0]].neighbors;
                for (int j = 0; j < (*pArbGraph).network[randObj[0]].neighbors; ++j) {
                    go += receive*((*pArbGraph).network[randObj[0]].neighbor[j]->attribute == 'c');
                }
                gn = - cost * ((*pArbGraph).network[randObj[0]].neighbor[randObj[1]]->attribute == 'c') * (*pArbGraph).network[randObj[0]].neighbor[randObj[1]]->neighbors;
                for (int j = 0; j < (*pArbGraph).network[randObj[0]].neighbor[randObj[1]]->neighbors; ++j) {
                    gn += receive*((*pArbGraph).network[randObj[0]].neighbor[randObj[1]]->neighbor[j]->attribute == 'c');
                }
                double possibility = sigmoid((gn - go)/.1);
                if (1.*abs(rand())/RAND_MAX < possibility) {
                    (*pArbGraph).network[randObj[0]].attribute = (*pArbGraph).network[randObj[0]].neighbor[randObj[1]]->attribute;
                }
            }
        }
        density[0] = density[1];
        for (int i = 0; i < (*pArbGraph).nodes; ++i) {
            density[1] += 1.*(((*pArbGraph).network)[i].attribute == 'c');
        }
        density[1] /= 1.*(*pArbGraph).nodes;
        iterCount++;
    } while (fabs(density[0] - density[1]) > 1.e-6 && iterCount < (*pArbGraph).nodes);
    return (density[1]);
}

double mcIntegral(int mcSteps, double start, double end, double* constantParameters, double (*pFunction)()) {
    double igv = 0.;
    for (int i = 0; i < mcSteps; ++i) {
        double loc = (start < end ? start : end) + fabs(end - start)*abs(rand())/RAND_MAX;
        igv += pFunction(loc, constantParameters)/mcSteps;
    }
    igv *= fabs(end - start);
    return igv;
}
